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I.  INTRODUCTION 


Over  the  past  several  years.  Mission  Research  Corporation  (MRC)  has 
developed  two  sets  of  general-purpose  computer  programs  for  multidimensional, 
electromagnetic  simulations  of  relativistic  beams  and  plasmas.  Applications 
have  included  magnetic  insulation,  electron  diodes,  electron  and  ion  acceler¬ 
ation,  electron-layer  magnetic  compression,  ion-focused  transport  of  electron 
beams,  beam-plasma  instabilities,  electron  beam  atmospheric  propagation,  and 
microwave  generation.  In  addition,  MRC  has  several  simpler,  faster  running 
codes  for  investigating  electron  beam ^equi 1 ibrium  and  linear  stability  proper¬ 
ties.  For  the  most  part,  these  smaller  codes  have  been  employed  in  acceler¬ 
ator  research.  The  simulation  codes,  together  with  the  more  important 
equilibrium  and  stability  codes,  are  listed  in  Table  1. 

This  report  summarizes  the  capabilities  and  describes  the  physical  models 
of  these  codes,  especially  as  they  bear  upon  microwave  source  investigations. 
Section  II  presents  briefly  the  principal  capabilities  of  each  of  the  fourteen 
codes  in  Table  1.  Section  III  begins  with  an  overview  of  particle-in-cel 1 
(PIC)  simulation  techniques  and  then  enumerates  the  unique  aspects  of  each  of 
the  general  purpose  PIC  simulation  codes.  The  models  employed  in  the 
remaining  computer  programs  are  discussed  in  Section  IV, 

All  programs  listed  in  Table  1  except  CPROP  run  on  CRAY-1  computers  under 
the  CTSS  operating  system.  As  noted  in  Section  II,  some  run  on  other  com¬ 
puters  as  well.  Copies  reside  at  the  Los  Alamos  National  Laboratory  (LANL) 
computer  center  and  elsewhere.  The  U.  S.  Government  either  owns  or  has  a 
royalty-free  license  to  use  any  of  the  programs. 
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TABLE  1.  MAJOR  MRC  COMPUTER  PROGRAMS  FOR  MICROWAVE 
AND  CHARGED  PARTICLE  BEAM  RESEARCH 


CODE 

APPLICATIONS 

REFERENCE 

CCUBE* 

2-D  PIC  Microwave  and  Beam  Simulations 

1 

IVORY 

3-D  PIC  Microwave  and  Beam  Simulations 

2 

ISIS* 

2-D  PIC  Microwave  and  Beam  Simulations 

3 

MAGIC 

2-D  PIC  Microwave  and  Beam  Simulations 

4 

SOS 

3-D  PIC  Microwave  and  Beam  Simulations 

5 

CPROP 

2-D  PIC  Beam  Atmospheric  Propagation 

6 

IPROP 

3-D  PIC  Beam  Atmospheric  and  IFR  Propagation 

7 

SCRIBE 

2-D  Electron  Diode  Equilibria 

8 

ARCTIC 

2-D  Ion-Focused  Transport  Equilibria 

9 

BALTIC 

Beam  Transverse  Stability  in  Accelerators 

10 

GRADR* 

Laminar  Beam  Stability,  Microwave  Growth 

11 

ORBIT 

1-D  Beam  Kinetic  Equilibria 

12 

KMRAD 

Kinetic  Beam  Stability,  Microwave  Growth 

13 

BTRSQ 

Recirculating  Accelerator  Instabilities 

14 

♦Developed  by  Los  Alamos  National  Laboratory  and  used  with  its 

permission. 
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II.  CODE  CAPABILITIES 


1.  OVERVIEW 

PIC  codes  treat  the  detailed  electromagnetic  interactions  among  one  or 
more  plasmas,  charged  particle  beams,  and  metallic  structures.  They  have  been 
employed  fruitfully  in  every  branch  of  plasma  physics  from  fusion  energy  to 
astrophysics;  the  linear  and  nonlinear  evolution  of  instabilities  is  a  common 
theme.  In  high-current  electron  accelerator  research,  PIC  codes  are  used  to 
treat  pul sed-power  flow  in  magnetically  insulated  transmission  lines,  electron 
beam  generation  in  complex  diodes,  beam  transport  and  inductive  acceleration 
in  vacuum  drift  tubes,  and  beam  transport  in  low  density  ion  channels.  Beam 
propagation  in  partially  ionized  gases  is  a  related  application.  High-power 
microwave  sources  are  modeled  effectively  by  PIC  codes,  which  can  follow 
simultaneously  the  electron  beam  nonlinear  dynamics  and  the  resonant  electro¬ 
magnetic  response  of  the  microwave  cavity.  Microwave  extraction  also  can  be 
examined.  The  microscopic  treatment  provided  by  PIC  codes  is  especially  valu¬ 
able  in  verifying  the  frequency  and  phase  coherence  of  extended  microwave 
sources.  PIC  codes  also  can  be  utilized  in  microwave  antenna  design,  propaga¬ 
tion  analysis,  and  coupling  studies,  although  simpler  codes  typically  are 
preferable. 

2.  CCUBE 

The  PIC  code  CCUBE  is  used  extensively  to  model  microwave  generation 
(Refs.  15-18)  and  beam  acceleration  (Refs.  19-22).  Although  only  two- 
dimensional,  the  code  is  very  flexible  in  that  it  can  run  in  any  orthogonal 
coordinate  system.  Injection  and  extraction  of  particles  and  fields  at 
boundaries  allows  it  to  treat  complicated,  open-ended  microwave  devices.  The 
code  has  extensive  diagnostics,  including  particle  and  field  plots,  time 
histories  of  various  quantities,  and  color  movies.  A  graphics  post-processor 
package,  PEGASUS,  is  available  (Ref.  23).  CCUBE  runs  on  VAX  780  and  CONVEX 
computers,  as  well  as  on  CRAY-1  computers. 
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3.  IVORY 


IVORY  is  a  three-dimensional  generalization  of  CCUBE.  Field  variation 
in  the  third  dimension  is  treated  by  Fourier  series  of  one  to  about  six  terms. 
As  a  consequence,  the  code  is  best  suited  for  problems  in  which  the  particle 
dynamics  are  only  weakly  nonlinear  in  the  third  dimension  and  boundary  condi¬ 
tions  are  approximately  two-dimensional.  Examples  include  free  electron 
lasers,  many  microwave  sources,  particle  beam  transport  (Refs.  24-27),  and 
plasmoid  propagation  across  a  magnetic  field  (Ref.  28).  IVORY  is  very 
efficient  for  such  problems,  requiring  computer  resources  only  a  few  times 
greater  than  those  needed  for  corresponding  two-dimensional  problems. 


4.  ISIS 


ISIS  is  a  significantly  enhanced,  but  still  two-dimensional,  version  of 
CCUBE.  Special  features  include  Monte  Carlo  transport  of  high  energy  par¬ 
ticles,  volume  creation  of  particles  by  high  energy  electrons  and  photons,  and 
a  hydrodynamic  treatment  of  imploding  resistive  liners.  The  layout  of  compli¬ 
cated  electrode  geometries  is  largely  automated,  and  a  new  algorithm  for 
improved  treatment  of  curved  metallic  surfaces  is  under  development.  It  and 
CCUBE  are  used  at  LANL  in  vircator  (Refs.  29-30),  electron  diode,  collective 
acceleration,  ring  compression,  and  beam-plasma  interaction  studies. 


5.  MAGIC 


MAGIC  is  a  two-dimensional  PIC  code  quite  similar  in  capabilities  to 
CCUBE  but  developed  independently.  It  is  employed  at  Sandia  National 
Laboratories  (SNL)  and  elsewhere  for  magnetic  insulation,  electron  diode 
(Ref.  31),  collective  acceleration  (Ref.  32),  beam  transport  (Ref.  33),  and 
microwave  generation  studies  (Ref.  34). 


6.  SOS 

SOS  is  a  fully  three-dimensional ,  general  purpose  simulation  code.  It 
accommodates  essentially  arbitrary  boundary  conditions  in  cartesian  or 
cylindrical  coordinates.  Wave  and  particle  emission  and  absorption,  metallic 
structures,  dielectrics,  and  air  conductivity  packages  are  included.  Sophis¬ 
ticated  data  management  techniques  allow  it  to  treat  problems  of  immense  size 
and  complexity,  limited  only  by  the  availability  of  computer  time.  A  compact, 
"user  friendly"  data  entry  system  simplifies  problem  specification.  In  addi¬ 
tion  to  plasma  physics  research,  SOS  often  is  used  in  nuclear  weapons  effects 
studies.  SOS  runs  on  VAX  780  and  CONVEX  computers,  as  well  as  on  CRAY-1 
computers. 

7.  CPROP 

CPROP  is  an  enhanced  and  optimized  version  of  CCUBE  for  axisymmetric 
investigations  of  electron  beam  propagation  in  air.  Special  features  include 
a  mesh  which  moves  with  the  beam,  a  p'rtially  implicit  field  solver  to  permit 
large  time  steps,  a  multispecies  air  chemistry  package,  and  a  delta-ray  crea¬ 
tion  and  destruction  procedure.  CPROP  is  not  limited  by  the  paraxial  particle 
and  frozen  field  approximations  used  in  many  other  propagation  codes.  The 
code  has  been  employed  in  Nordsieck  expansion  (Ref.  35),  beam  front  erosion, 
and  hollowing  instability  studies  supporting  RADLAC  (Refs.  36-37),  IBEX 
(Refs.  38-40),  PHERMEX  (Ref.  41),  and  FX-100  (Ref.  42)  experiments.  CPROP 
runs  on  CDC-7600  computers  under  the  LTSS  operating  system. 

8.  IPROP 


IPROP  is  a  three-dimensional  electron  beam  atmospheric  propagation  code 
obtained  by  adding  the  special  features  of  CPROP  to  IVORY.  Diamagnetic 
effects  are  included  to  treat  the  RADLAC  beam  rotation  properly.  The  code  can 
model  all  resistive  instabilities  --  hollowing,  hose,  and  f i 1 amentati on  — 
either  separately  or  in  combination  (Ref.  36).  Two  major  channel-tracking 
discoveries  were  made  with  IPROP  in  1986  (Refs.  43,  44).  With  the  air 
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conductivity  package  omitted,  I  PROP  also  is  used  for  research  on  electron  beam 
transport  in  low-density  ionized  channels,  the  so-called  ion-focused  regime 
(IFR)  (Ref.  45).  Research  centers  on  beam  front  erosion,  two-stream 
instabilities,  and  channel  electron  trapping. 

9.  SCRIBE 

SCRIBE  can  be  described  as  a  two-dimensional,  steady-state  PIC  code.  It 
is  based  on  Hermannsfeldt ‘ s  Electron  Trajectory  Program  (Ref.  8).  It  computes 
axisymmetric  electron  flow  patterns  in  sel f-consi stently  determined  static 
electric  and  magnetic  fields.  Being  time-independent,  the  code  is  fast  and 
relatively  free  of  numerical  noise.  SCRIBE  is  ideal  for  studying  electron 
equilibrium  behavior  in  diodes  (Ref.  46),  accelerators  and  microwave  tubes.  A 
particularly  noteworthy  achievement  with  SCRIBE  is  the  design  of  very  low 
emittance  electron  diodes  and  transport  systems  for  free  electron  laser  stud¬ 
ies  at  the  Naval  Research  Laboratory  (Ref.  47).  SCRIBE  runs  on  VAX  780 
computers,  as  well  as  on  CRAY-1  computers. 

10.  ARCTIC 

ARCTIC  is  similar  to  SCRIBE  but  runs  in  the  beam  frame  and  solves  the 
frozen  field  equations.  It  was  written  to  determine  electron  and  ion  trajec¬ 
tories  and  associated  steady-state  fields  for  relativistic  electron  beam 
propagation  in  IFR  channels.  For  such  applications  ARCTIC  is  far  more  effi¬ 
cient  than  alternative  codes,  especially  when  very  long  beam  pulses  and  large 
drift  tubes  are  involved. 

11.  BALTIC 

BALTIC  is  a  simple  code  used  to  follow  the  transverse  oscillations  of 
the  centroid  of  a  particle  beam  in  an  accelerator.  Specifically,  it  computes 
the  linear  growth  of  the  resistive  wall  (Ref.  48),  image  displacement  (Ref. 
24),  and  beam  breakup  instabilities  (Refs.  10,  49)  for  long  pulse  beams  in 
multigap  accelerators  like  RADLAC.  With  minor  modifications  it  could  be 
employed  to  predict  gain  in  Scantron  devices. 


12.  GRADR 


GRADR  is  a  general-purpose  equilibrium  and  linear  dispersion  code  for 
laminar,  cyl indrical ly  symmetric,  radially  inhomogeneous  particle  beams.  Con¬ 
ducting  boundaries  and  resistive  and  dielectric  mediums  are  allowed.  Starting 
from  user-supplied  constraint  equations,  the  code  first  computes  the  beam 
equilibrium  current  and  field  profiles.  It  then  determines  wave  frequencies 
and  growth  rates  for  arbitrary,  user-requested  wave  numbers.  GRADR  is  inter¬ 
active  and  exceedingly  fast.  Typically,  the  code  is  used  in  collective  accel¬ 
eration  (Refs.  19,  50-52)  and  beam  stability  studies  (Refs.  48,  53-56).  With 
modifications,  it  could  be  employed  for  microwave  generation  in  slow  wave 
structures  (Ref.  57). 

13.  ORBIT 

ORBIT  and  KMRAD  are  used  in  place  of  GRADR  for  nonlaminar  beams,  the  more 
common  situation.  ORBIT  determines  the  equilibrium  radial  profiles  of  the 
particle  currents  and  fields,  based  on  a  user-provided  subroutine  specifying 
the  particle  distribution  function  in  terms  of  constants  of  the  motion.  A 
family  of  distribution  function  subroutines  for  relativistic  Bennett  beams  and 
for  radial  diodes  is  available  (Ref.  58). 

14.  KMRAD 

KMRAD  is  a  three-dimensional  linearized  PIC  code.  A  cyl indrical ly 
symmetric,  radially  inhomogeneous  particle  distribution  is  specified  by  the 
user,  possibly  with  the  aid  of  ORBIT.  (An  interface  between  these  codes  is 
provided.)  It  then  computes  the  temporal  evolution  of  perturbed  particle 
orbits  and  electromagnetic  fields  for  arbitrary  wave  numbers  chosen  by  the 
user.  The  code  is  interactive  and  moderately  fast.  KMRAD  is  used  regularly 
for  parametric  studies  of  instability  growth  for  propagating  particle  beams  in 
air  (Refs.  36,  59-61)  and  in  IFR  channels.  Smooth-bore  magnetron,  Cherenkov 
masers,  gyrotrons,  and  orbitrons  also  can  be  studied.  With  modifications,  the 
code  can  be  employed  for  microwave  generation  in  slow  wave  structures. 
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15.  BTRSQ 
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BTRSQ  is  the  most  detailed  of  a  family  of  dispersion-relation  solvers  for 
high-current  betatrons  (Refs.  25,  62).  The  beam  is  treated  in  the  rigid  disk 
limit,  while  the  fields  are  determined  exactly  for  a  cavity  of  rectangular 
minor  cross  section.  The  code  is  interactive  and  exceedingly  fast.  Many  new 
properties  of  the  negative  mass  instability  in  modified  betatrons  have  been 
discovered  recently  with  BTRSQ  (Refs.  14,  27,  63). 


Oispersion  relation  codes  for  microwave  generation  can  be  developed 
rapidly  as  the  needs  arise.  A  convenient  shell,  developed  originally  for 
GRADR,  provides  user  interaction,  error  recovery,  a  Muller's  method  (Ref.  64) 
complex  root  evaluater,  and  a  graphics  interface.  The  programmer  merely 
writes  a  dispersion  function  evaluation  routine  for  the  new  device.  For 
instance,  the  Pierce  diode  electromagnetic  stability  code  PEM2D  being  written 
for  AFWl  by  MRC  uses  this  shell. 


III.  PIC  CODE  ALGORITHMS 


1.  GENERAL  FEATURES  OF  PIC  SIMULATION  CODES 

PIC  codes  determine  the  time  evolution  of  complex  plasmas  by  computing 
the  dynamics  of  many  thousands  of  representative  plasma  particles  (electrons 
and/or  ions)  moving  in  electromagnetic  fields  externally  applied  or  produced 
by  the  plasma  itself.  Thus,  PIC  codes  provide  the  most  fundamental  and 
detailed  representation  possible  of  plasma  problems.  In  effect,  they  solve 
the  Vlasov  equation.  Of  course,  this  precision  comes  at  the  cost  of  substan¬ 
tial  computer  requirements,  and  for  this  reason  PIC  codes  should  be  employed 
only  when  simpler  numerical  or  analytical  techniques  are  inadequate. 

The  electromagnetic  fields  are  defined  on  a  regular  mesh  in  one  two,  or 
three  dimensions,  depending  on  the  symmetry  of  the  problem  to  be  solved.  The 
mesh  can  be  in  rectangular,  cylindrical,  or  other  desired  geometry.  At  each 
time  step,  new  electric  (E)  and  magnetic  (B)  fields  are  computed  by  advancing 
the  finite  difference  approximations  to  Maxwell's  equations, 


it  v  x  B  -  J 


=  -V  x  E  (2) 

using  currents  (J)  determined  from  the  plasma  particle  motion  on  the  previous 
time  step.  Alternati vely,  equations  for  the  scaler  and  vector  potentials  can 
be  solved.  Boundary  conditions  are  required  to  define  spatial  derivatives  at 
the  mesh  edges.  Wave  reflecting  (i.e.,  metallic)  or  periodic  boundary  condi¬ 
tions  are  common  choices.  More  complicated  boundaries  allow  electromagnetic 
waves  to  be  launched  into  the  computational  region  or  to  leave  it.  Metal 
structures  at  locations  within  the  computational  region  are  represented  simply 
by  setting  electric  field  components  to  zero  there.  Representing  dielectrics, 
resistors,  and  simple  driven  antennas  is  only  slightly  more  difficult. 
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Particle  momenta  (P  =  yV),  positions  (X),  and  energy  ( y)  are  then 
advanced  using  the  relativistic  equations  of  motion  with  the  newly  computed 
fields. 

f  sE,'-»  (3) 

Hr  *v  <4> 

Y  =  (P2  +  1)1/2  (5) 

(Replace  y  by  1  for  nonrelati vistic  problems.)  The  fields  appearing  in  Equa¬ 
tion  3  are  those  at  the  particle  location,  obtained  from  the  fields  at  nearby 
mesh  points  by  (typically,  linear)  interpolation.  When  a  particle  leaves  the 
computational  region,  it  is  destroyed  or  it  is  returned  to  the  mesh  by  some 
prescribed  procedure  (e.g.,  reflected).  By  the  same  token,  particles  can  be 
injected  from  boundaries,  a  feature  particularly  useful  in  particle  beam  simu¬ 
lations.  After  a  particle's  new  position  and  momenta  have  been  determined, 
its  contribution  to  the  plasma  currents  is  obtained  by  interpolating  V  =  P/y 
to  nearby  mesh  points. 

This  cycle  of  advancing  fields  based  on  particle  currents  and  then 
advancing  particles  based  on  the  new  fields  is  repeated  hundreds  of  times  in  a 
typical  simulation.  The  time  step  is  set  by  the  smallest  time  scale  in  the 
calculation,  which  may  be  the  plasma  oscillation  period,  the  electron  cyclo¬ 
tron  period,  or  the  Courant  time  (of  order  the  time  for  a  light  wave  to  cross 
a  cell  in  the  mesh).  Progress  has  been  made  in  the  last  few  years  at  sur¬ 
mounting  the  Courant  limit,  which  is  numerical  rather  than  physical  in  charac¬ 
ter  (Ref.  65).  Cell  dimensions  must  be  small  compared  to  spatial  scales  of 
interest. 

The  PIC  code  running  times  and  memory  requirements  are  highly  problem- 
dependent.  The  CPU  times  on  a  CRAY-1  computer  typically  range  between  15  min 
and  4  h,  although  20-h  runs  are  not  unheard  of.  Corresponding  central  memory 
needs  vary  between  2*10^  and  4 •  106  words.  At  least  two  fast,  large  capacity 
disks  also  are  needed.  Historically,  the  physics  problems  attempted  with  PIC 
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codes  have  expanded  to  consume  the  maximum  resources  available  in  each 
generation  of  computers. 


PIC  codes  usually  have  extensive  graphics  output  capabilities  and 
operating-system  interfaces.  Advances  in  PIC  technique  are  described  in 
several  books  (Refs.  66-67). 

2.  CCUBE 

CCUBE  (Ref.  1)  is  a  general-purpose  two-dimensional  PIC  code,  and  as  such 
operates  along  the  lines  described  above.  It  does,  however,  have  a  number  of 
special  features  which  increase  its  flexibility.  The  Galerkin  finite-element 
algorithm  for  advancing  particle  quantities  (Ref.  68)  eliminates  the  numerical 
Cherenkov  instability  (Ref.  69)  and  ameliorates  some  other  numerical  prob¬ 
lems.  The  electromagnetic  fields  incorporate  backward  biasing  to  damp 
unwanted  high-frequency  field  fluctuations  and  to  further  suppress  possible 
numerical  instabilities  (Ref.  70).  The  MRC-developed  time-biased  field  solver 
also  can  significantly  relax  the  Courant  time  step  constraint,  reducing  com¬ 
puting  costs.  Marder's  procedure  for  reducing  charge  conservation  errors  is 
being  added  (Ref.  71).  CCUBE  is  written  for  arbitrary  orthogonal  coordinates. 
The  user  can  convert  from,  say,  cylindrical  to  toroidal  coordinates  merely  by 
revising  routines  which  define  metric  elements.  Nonuniform  zoning  of  the 
spatial  mesh  is  accomplished  with  equal  ease. 

A  number  of  other  features  deserve  mention.  For  compatibility  with 
nonuniform  zoning,  variable  particle  weighting  within  particle  type,  with 
several  particle  species  permitted,  is  implemented.  The  code  supports  peri¬ 
odic,  wave  transmitting,  and  inhomogeneous  Dirichlet-Neumann  field  boundary 
conditions.  Particles  can  be  absorbed,  reflected,  or  injected  from  surfaces 
not  necessarily  coinciding  with  edges  of  the  mesh.  The  physics  for  wave 
launching  and  particle  field-emission  from  surfaces,  which  is  necessary  to 
accurately  model  microwave  and  particle  beam  sources,  is  incorporated.  Graph¬ 
ics  output  can  be  produced  either  directly  from  CCUBE  or  interactively  from 
the  PEGASUS  (Ref.  23)  postprocessor.  Options  include  microfiche  plots  of 
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various  slices  through  particle  phase  space,  particle  distribution  functions, 
contour  plots  and  one-dimensional  cuts  of  fields  and  currents,  and  histories 
of  particle  and  field  energies  and  other  selected  quantities.  Color  movies  of 
particle  and  field  data  can  be  generated  concurrently.  In  addition  to  the 
more  traditional  simulation  diagnostics  just  mentioned,  CCUBE  contains  numer¬ 
ical  equivalents  of  such  experimental  diagnostics  as  Faraday  cups,  calorim¬ 
eters,  compensated  diamagnetic  loops,  Rogowski  coils,  and  local  probes  for 
field  and  current  measurements.  Each  is  time  dependent  and  can  be  Fourier 
analyzed. 

3.  IVORY 

IVORY  (Ref.  2)  is  a  three-dimensional  generalization  of  CCUBE  and  shares 
with  it  most  of  the  special  features  described  above.  For  economy  of  opera¬ 
tion  and  clarity  of  interpretation,  variations  of  field  quantities  in  the 
third  (usually,  azimuthal)  coordinate  are  represented  by  spectral  rather  than 
finite  difference  methods;  i.e.,  the  fields  are  Fourier  decomposed  in  that 
dimension.  In  this  way  the  behavior  of  selected  modes  can  be  examined  without 
wasting  time  and  storage  on  other,  uninteresting  modes.  It  should  be  empha¬ 
sized  that  the  spectral  method  is  not  a  linearization  in  either  particle  or 
field  behavior.  Of  course,  the  greater  the  nonlinearity  of  a  problem,  the 
more  azimuthal  modes  are  needed  to  resolve  it  accurately. 

Thus,  IVORY  is  particularly  suitable  for  configurations  which  depart 
only  weakly  from  symmetry  in  one  coordinate,  such  as  a  particle  beam  in  an 
accelerator  or  microwave  device.  The  code  provides  three-dimensional  results 
at  costs  comparable  to  two-dimensional  simulations  in  such  cases. 

4.  ISIS 

The  two-dimensional  ISIS  code  (Ref.  3)  also  is  a  derivative  of  CCUBE, 
modified  to  allow  more  accurate  treatment  of  particle  beam  interaction  with 
materials.  Electron-  and  photon-generated  delta  rays  can  be  created,  trans¬ 
ported,  and  destroyed  in  the  simulations.  A  Monte  Carlo  treatment  of  Moliere 
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scattering  is  provided  for  particles  passing  through  foils  and  other 
materials.  Energy  loss  is  modeled  by  Bethe's  continuous-sl owing-down  approxi¬ 
mation.  A  resistive  magnetohydrodynamic  representation  of  metallic  shells  is 
useful  in  magnetic  compression  studies.  Bulk  resistivity  terms  have  been 
added  to  the  field  equations. 

Several  numerical  enhancements  are  noteworthy.  Bilinear  interpolation  of 
currents  and  fields  is  permitted  in  addition  to  the  Galerkin  scheme.  A  cur¬ 
rent  lay-down  procedure  explicitly  guaranteeing  charge  conservation  is 
provided.  The  backward-biased  field  algorithm  is  inverted  exactly  rather  than 
iteratively.  The  algorithm  for  space-charge-limited  particle  emission  from 
electrodes  is  more  accurate,  and  an  improved  treatment  of  curved  electrodes  is 
being  developed.  The  procedure  for  specifying  the  locations  of  complex 
electrodes  on  the  mesh  is  automated. 

5.  MAGIC 

The  two-dimensional  code  MAGIC  (Ref.  4)  was  developed  almost  entirely 
independently  of  CCUBE.  Nonetheless,  the  two  are  very  similar  not  only  in 
their  features  but  also  in  their  internal  algorithms. 

Differences  are  noted  here.  Field  interpolation  in  MAGIC  is  bilinear, 
and  current  interpolation  is  explicitly  charge-conserving.  Triangular  zoning 
at  conductor  surfaces  improves  field  accuracy  there.  Wave  absorption  at 
boundaries  by  graded  electric  and  magnetic  resistivity  is  accommodated.  Input 
data  processing  is  a  bit  more  flexible.  There  is  no  postprocessor. 

6.  SOS 

SOS  (Ref.  5)  is  a  general  purpose  three-dimensional  PIC  code.  Unlike 
IVORY,  it  treats  field  variations  in  all  three  coordinates  by  finite  differ¬ 
ences.  The  code  was  originally  developed  for  system-generated  electromagnetic 
pulse  computations  but  has  been  generalized  to  treat  particle  beam,  microwave, 
and  magnetic  insulation  problems. 


A  distinguishing  feature  of  SOS  is  its  data  management  procedure.  At 
each  time  step,  the  code  processes  the  three-dimensional  mesh  and  the  parti¬ 
cles  on  it  as  a  sequence  of  two-dimensional  slices.  Only  a  few  contiguous 
slices  are  in  computer  memory  at  any  instant,  while  the  rest  are  stored  on 
disk.  (In  most  other  PIC  codes,  fields  and  currents  are  retained  entirely  in 
computer  memory,  while  particles  are  cycled  between  disk  and  memory.)  In  this 
way,  extremely  large  computations  can  be  performed,  limited  only  by  budget. 

SOS  has  an  X-ray-generated  Compton  electron  creation  routine. 

Conductivity  generation  by  Compton  electrons  and  other  particles  is  determined 
using  a  three-species  lumped  air  chemistry  model  with  direct  ionization,  ava¬ 
lanche,  and  recombination.  Bulk  Compton  currents  and  conductivity  in  dielec¬ 
trics  are  treated  by  a  special  model.  Exponential  integrating  factors  in  the 
field  algorithm  stably  accommodate  arbitrarily  variable  conductivity.  A 
family  of  subgrid  models  is  provided  for  structures  too  small  to  be  repre¬ 
sented  accurately  by  standard  finite  differences.  Other  features  of  the  codes 
are  similar  to  those  in  MAGIC. 

7.  CPROP 

CPROP  (Ref.  6)  is  a  two-dimensional  PIC  code  designed  for  investigating 
the  axisymmetric  dynamics  of  high-current  relativistic  electron  beam  propaga¬ 
tion  in  air.  It  is  based  on  CCUBE  and  thus  possesses  most  of  the  diagnostic 
and  other  capabilities  of  that  program.  Three  major  enhancements  were 
required  to  treat  beam  transport  in  air:  an  electromagnetic  field  solver 
accommodating  arbitrary  scalar  conductivity  and  a  moving  coordinate  mesh,  a 
conductivity  generation  package,  and  a  Moliere  electron  scattering  routine. 
CPROP  also  contains  a  high  energy  delta-ray  generation  and  transport  package. 

The  CPROP  field  algorithm  is  unique  in  several  respects  (Ref.  40). 
Calculations  are  carried  out  in  a  frame  moving  axially  at  a  user-specified, 
possibly  time-varying  velocity.  The  field  evaluation  is  partially  imolicit, 
so  that  the  Courant  limit  is  set  only  by  the  axial  cell  size.  The  so-called 
frozen  field  approximation  (Refs.  43,  72)  to  Maxwell's  equations  can  be 


treated  as  a  special  case,  if  desired.  Exponential  integrating  factors, 
similar  to  those  in  SOS,  accommodate  arbitrary  scalar  conducti vities.  Zoning 
is  nonuniform  both  axially  and  radially. 

Both  the  BMCOND  (Ref.  73)  and  the  PHOENIX  (Ref.  74)  air  chemistry  models 
are  available.  The  former  solves  rate  equations  for  gas  electron  density  and 
temperature  and  for  the  densities  of  several  ionic  species.  The  latter  is  a 
lumped  parameter  model  similar  to  that  in  SOS.  Moliere  scattering  (Refs.  35, 
75)  is  implemented  in  CPROP  by  applying  a  deflection  to  the  beam  electrons 
every  few  time  steps.  The  deflection  angle  is  chosen  randomly  from  a  set  of 
previously  computed  small  angles  forming  a  truncated  Moliere  distribution  or 
from  an  analytical  expression  for  occasional  large  angles. 

8.  IPROP 

IPROP  (Ref.  7)  is  a  version  of  IVORY  into  which  most  of  the  special 
features  of  CPROP  have  been  transplanted.  A  new  lumped  parameter  conductivity 
model  calibrated  against  the  SAIC  HICHEM  code  is  used  (Ref.  76).  Features  of 
the  electromagnetic  field  solver  and  the  conductivity  package  were  essential 
to  recent  MRC  channel  tracking  discoveries  (Refs.  43-44). 


IV.  ALGORITHMS  OF  OTHER  MRC  CODES 


The  remaining  computer  programs  are  much  faster  than  PIC  codes.  SCRIBE, 

ARCTIC,  8ALTIC,  and  KMRAD  typically  use  no  more  than  a  few  minutes  of  CRAY-1 

computer  time;  GRADR,  0R8IT,  and  BTRSQ  only  seconds.  None  of  these  codes 
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needs  more  than  10  words  of  memory. 

1.  SCRIBE 

When  only  a  steady-state  simulation  of  electron  beam  generation  and 
transport  is  needed,  great  economy  can  be  achieved  with  SCRIBE.  It  is  a  ver¬ 
sion  of  the  Stanford  Linear  Accelerator  (SLAC)  code  EGUN  (Ref.  8),  extensively 
modified  by  MRC  for  research  on  intense  beam  diodes.  SCRIBE  self-consistently 
solves  two-dimensional  problems  (usually  in  (r,z)  geometry)  by  an  iterative 
method.  First,  a  solution  to  Laplace's  equation  is  obtained,  then  current 
emission  is  calculated,  and  the  electron  trajectories  are  integrated  using 
the  fields.  The  trajectories  deposit  space  charge  at  the  mesh  points,  and  a 
solution  to  Poisson's  equation  is  calculated.  Current  emssion  is  then  recal¬ 
culated  and  the  cycle  repeated  until  the  potentials  and  trajectories  are 
unchanged  from  cycle  to  cycle.  Improved  algorithms  assure  rapid  convergence. 
All  relevant  self  and  applied  electric  and  magnetic  fields  (except  the  beam- 
generated  axial  magnetic  field)  are  included  in  the  calculations.  Both 
Dirichlet  and  Neumann  boundaries  are  available,  and  surfaces  are  not 
restricted  to  lie  along  the  grids. 

2.  ARCTIC 

ARCTIC  is  a  newly  developed  two-dimensional  code  for  determining  the 
steady-state  flow  of  channel  electrons  and  ions  for  IFR  beam  transport  prob¬ 
lems.  At  present,  the  beam  is  treated  as  a  rigid  rod  propagating  at  approxi¬ 
mately  the  speed  of  light.  Calculations  are  performed  in  the  beam  frame. 
Representative  channel  electron  and  ion  trajectories  are  computed  based  on 
self-consistent  fields,  and  the  corresponding  currents  interpolated  onto  a 
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spatial  mesh.  Because  fields  are  evaluated  in  the  frozen  field  limit  (Refs. 
43,  72),  no  iterations  are  required  to  reach  the  solution.  In  many  respects, 
ARCTIC  resembles  a  one-dimensional  PIC  code  in  which  the  axial  coordinate 
takes  the  place  of  time. 

A  three-dimensional  version  of  ARCTIC  is  planned  for  the  future.  Repre¬ 
senting  the  beam  cross  section  by  an  envelope  model  also  is  attractive, 
although  iterative  solution  of  the  resulting  system  of  equations  would  then  be 
necessary. 

3.  BALTIC 

BALTIC  (Ref.  10)  provides  a  less  detailed  but  immensely  faster  treatment 
than  IVORY  of  transverse  beam  instabilities  in  accelerators.  Presently,  it 
evaluates  the  simultaneous  evolution  of  resistive  wall  (Ref.  48),  image  dis¬ 
placement  (Ref.  24),  and  beam  breakup  (Refs.  10,  49)  instabilities  in 
multiple-gap  linear  induction  and  radio  frequency  accelerators.  Arbitrary 
spatial  variations  in  the  beam  energy  and  axial  magnetic  field  strength  are 
allowed.  Provision  exists  for  external  forces.  Adding  centrifugal  force  and 
a  vertical  magnetic  field  for  recirculating  accelerator  studies  is 
straightforward. 

BALTIC  employs  a  number  of  simplifying  assumptions  for  computational 
efficiency.  The  beam  is  represented  as  a  chain  of  rigid  disks  propagating 
forward  at  a  uniform  speed.  Electric  and  magnetic  fields  entering  into  the 
equations  of  transverse  motion  are  computed  in  the  long  wavelength  limit,  and 
accelerating  gaps  are  treated  in  the  thin  lens  approximation.  Coupling  coef¬ 
ficients  entering  into  various  instabilities  can  be  estimated  analytically 
(Ref.  77)  or  determined  from  experiments  or  more  detailed  numerical 
calculations. 


4. 


GRADR 


GRADR  (Ref.  11)  determines  the  equilibrium  profiles  and  linear  eigenmodes 
of  cyl indrical ly  symmetric,  radially  inhomogeneous,  relativistic  electron 
beams  in  the  laminar  flow  approximation.  The  cylindrical  equilibrium  is  eval¬ 
uated  from  six  nonlinear,  coupled,  algebraic  and  first-order  differential 
equations.  Four  follow  from  Maxwell's  equations  plus  the  fluid  equations, 
while  two  can  be  specified  arbitrarily  to  select  a  desired  equilibrium  (Ref. 
78).  The  two  constraint  equations  should,  of  course,  be  physically  realiz¬ 
able.  One  constraint  almost  always  is  appropriate:  total  particle  energy, 
kinetic  plus  potential,  must  be  constant  across  the  beam  and  equal  to  the 
injection  energy.  As  a  second  constraint,  typically  we  let  the  current  den¬ 
sity  profile  of  the  beam  be  specified  at  injection.  Conservation  of  canonical 
angular  momentum  then  determines  the  current  density  profile  within  the  accel¬ 
erator  drift  tube.  The  resulting  system  of  equations  is  solved  iteratively. 

Given  an  equilibrium  determined  in  this  or  any  other  way,  GRADR  then 
solves  the  corresponding  linearized  (small  amplitude  wave)  equations  to  obtain 
eigenmodes  and  eigenvalues.  The  linear  equations  form  a  fourth-order  differ¬ 
ential  system  in  radius  and  are  integrated  by  a  Fehlberg  fourth-fifth-order 
Runge-Kutta  routine  (Ref.  79).  Muller's  method  is  employed  in  finding  wave 
frequencies  and  growth  rates  (Ref.  64).  The  options  of  computing  wave  energy 
and  the  adiabatic  variation  of  wave  amplitude  with  changing  beam  parameters 
are  available  (Ref.  51).  Including  background  plasma  conductivity  in  GRADR 
allows  study  of  laminar  beam  resistive  instabilities  in  air  (Ref.  56).  Modi¬ 
fying  GRADR  boundary  conditions  to  accommodate  slow-wave  structures  is 
straightforward  (Ref.  57)  and  would  allow  gain  calculations  for  many  microwave 
devices. 

5.  ORBIT 

ORBIT  (Ref.  12)  calculates  one-dimensional,  radially  inhomogeneous, 
cyl indrical ly  symmetric  beam  and  plasma  equilibria.  The  code  solves  the 
relativistic  VI asov-Maxwel 1  equations  for  user-supplied  particle  distribution 
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functions.  Distributions  are  specified  in  terms  of  particle  constants  of 
motion:  canonical  axial  momentum,  canonical  angular  momentum,  and  energy. 

The  kinetic  framework  includes  finite  temperatures  in  a  natural  manner  and 
allows  investigation  of  phenomena  associated  with  the  detailed  momentum-space 
distributions  that  cannot  be  treated  using  a  macroscopic  fluid  description. 

ORBIT  has  been  used  primarily  for  investigating  equilibria  in 
relativistic  electron  beams  and  magnetically  insulated  diodes.  Relativistic 
electron  beam  equilibria  with  both  sharp  and  smooth  edges  have  been  examined, 
and  distribution  functions  for  hollow  beams  identified.  For  magnetically 
insulated  diodes,  both  axial  and  azimuthal  magnetic  fields  have  been  exam¬ 
ined.  Finite  temperature  effects,  gas  prefill,  and  electron  injection  with 
nonzero  energy  have  been  investigated.  The  use  of  ORBIT  has  eliminated  making 
many  of  the  assumptions  and  approximations  necessary  in  the  analytic  treatment 
of  magnetic  insulation.  A  catalog  of  distribution  functions  and  their 
properties  has  been  compiled  for  use  with  ORBIT  (Ref.  58). 

6.  KMRAD 

KMRAD  (Ref.  13)  is  a  linearized  but  fully  relativistic  and 
electromagnetic  linear  stability  code  for  arbitrary  cyl indrical ly  symmetric 
plasma  systems.  It  consists  of  two  parts:  a  one-dimensional ,  radially 
resolved,  nonlinear  PIC  code;  and  a  three-dimensional,  spectrally  resolved, 
linear  PIC  code.  The  one-dimensional  routines  determine  particle  orbits  and 
corresponding  fields  in  a  cylindrical ly  symmetric,  possibly  slowly  evolving, 
equilibrium.  It  can  be  initialized  from  ORBIT,  and  an  interface  exists 
between  the  codes. 

The  three-dimensional  routines  compute  the  temporal  growth  of  small 
amplitude  perturbations  of  the  particle  equilibrium  orbits,  and  the  fields 
arising  from  these  perturbations.  The  three-dimensional  fields  are  Fourier 
decomposed  axially  and,  azimuthal ly,  and  one  mode  at  a  time  is  analyzed.  The 
instability  growth  rate  of  the  fastest  unstable  wave  of  a  selected  mode-number 
pair  (kz,m)  can  be  determined  from  the  exponential  growth  of  field  energy  in 


the  simulation.  User-specified  radially  inhomogeneous  background  conductivity 
and  return  current  profiles  allow  KMRAD  to  treat  resistive  instabilities.  The 
code  is  interactive  and  quite  fast.  In  effect,  it  simulates  three-dimensional 
linear  dynamics  but  requires  computer  resources  comparable  to  those  of  a  one¬ 
dimensional  code.  KMRAD  plus  ORBIT  constitutes  a  kinetic  beam  generalization 
of  GRADR.  Like  GRADR,  it  can  be  modified  to  handle  microwave  growth  in  slow 
wave  structures. 

7.  BTRSQ 

BTRSQ  (Ref.  14)  calculates  the  normal  mode  frequencies  and  growth  rates, 
if  any,  of  a  particle  beam  in  a  modified  betatron  or  similar  device.  It  does 
so  by  numerically  evaluating  the  roots  of  an  analytical  dispersion  relation. 
The  dispersion  relation  was  derived  by  treating  the  beam  as  a  string  of  rigid 
disks  traveling  around  the  toroidal  acceleration  cavity  while  executing  small 
toroidal  and  poloidal  oscillations.  The  electromagnetic  fields  which  interact 
with  these  oscillations  were  evaluated  exactly  for  a  cavity  of  rectangular 
minor  cross  section  by  means  of  Green's  functions. 

Many  features  of  BTRSQ  are  patterned  after  those  in  GRADR.  It  has 
interactive  NAMELIST  input,  an  optimized  Muller's  method  root  finder,  error 
recovery  and  user-interrupt  procedures,  and  graphical  output.  It  is  extremely 
fast  and  very  easy  to  use.  Other  linear  stability  problems  can  be  addressed 
in  the  same  way  simply  by  providing  the  appropriate  dispersion  function 
routines.  This  has  been  done  in  several  instances. 
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